/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		displacement.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-04 17:48
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:26
*
*   Description:
*
=============================================================================*/
#ifndef DISPLACEMENT_H
#define DISPLACEMENT_H

#include "../../geometry/shape/nodal_dof.hpp"
#include "../../model/solve/assemble_elemental_mat.hpp"

namespace geoxfem {

namespace model {
struct displacement {};

template <>
template <typename DimensionType, typename degree>
MatrixXR assemble_gauss_kay<displacement, displacement> /*{{{*/
  ::apply(const geometry::ShapeInterface<DimensionType, degree>& shp,
          const GaussPoint<DimensionType>&                       gp) {
  MatrixXR dee = DeeMatrix(gp);
  MatrixXR bee = BeeMatrix(shp, gp);
  MatrixXR kay = bee.transpose() * dee * bee;
  GEOXFEM_ASSERT(!kay.hasNaN());
  return kay;
} /*}}}*/

template <>
template <typename DimensionType, typename degree>
VectorXR assemble_sub_rhs_for_tstep<displacement, displacement> /*{{{*/
  ::apply(const geometry::ShapeInterface<DimensionType, degree>& shp) {
  using Shape = const geometry::ShapeInterface<DimensionType, degree>;

  return algorithm::integral(shp, [](const Shape& shp1, const auto& gp) {
    MatrixXR bee    = BeeMatrix(shp1, gp);
    VectorXR sigma  = gp.getSigma().Voigt();
    VectorXR e_load = bee.transpose() * sigma;
    GEOXFEM_ASSERT(!e_load.hasNaN());
    return e_load;
  });
} /*}}}*/

template <>
template <typename DimensionType, typename degree>
VectorXR assemble_sub_rhs_for_iter<displacement, displacement> /*{{{*/
  ::apply(const geometry::ShapeInterface<DimensionType, degree>& shp) {
  using Shape = const geometry::ShapeInterface<DimensionType, degree>;

  return algorithm::integral(shp, [](const Shape& shp1, const auto& gp) {
    MatrixXR bee    = BeeMatrix(shp1, gp);
    VectorXR sigma  = gp.getSigmaAtTimeStep().Voigt();
    VectorXR e_load = bee.transpose() * sigma;
    GEOXFEM_ASSERT(!e_load.hasNaN());
    return e_load;
  });
} /*}}}*/

}  // namespace model

namespace geometry {

template <>
template <typename Vertex>
void set_nodal_dofs<model::displacement>::apply(Vertex& v) { /*{{{*/
  constexpr size_t dimension = traits::dimension<Vertex>::value;
  size_t           node      = v.getIndex();
  for(size_t i = 0; i != dimension; ++i)
    v.append(std::string("x") + std::to_string(i), model::Dof{node});
  if(v.ifXfem())
    for(size_t i = 0; i != dimension; ++i)
      v.append(std::string("Mx") + std::to_string(i), model::Dof{node});
  return;
} /*}}}*/

template <>
template <typename Vertex>
std::vector<model::Dof*> get_nodal_dofs<model::displacement> /*{{{*/
  ::apply(Vertex& v) {
  constexpr size_t         dimension = Vertex::dimension;
  std::vector<model::Dof*> dof_list;
  for(size_t i = 0; i != dimension; ++i)
    dof_list.emplace_back(v.getDof(std::string("x") + std::to_string(i)));
  if(v.ifXfem())
    for(size_t i = 0; i != dimension; ++i)
      dof_list.emplace_back(v.getDof(std::string("Mx") + std::to_string(i)));
  return dof_list;
} /*}}}*/

template <>
template <typename DimensionType, typename degree>
std::vector<model::Dof*> get_elemental_dof<model::displacement> /*{{{*/
  ::apply(const geometry::ShapeInterface<DimensionType, degree>& shp) {
  constexpr size_t         dimension = traits::dimension<DimensionType>::value;
  std::vector<model::Dof*> doflist;
  for(auto& v : neighborsIn(shp, 0))
    for(size_t i = 0; i != dimension; ++i) {
      std::string dofname = std::string("x") + std::to_string(i);
      doflist.emplace_back(v->toVertex().getDof(dofname));
    }
  if(shp.ifXfem())
    for(auto& v : neighborsIn(shp, 0))
      for(size_t i = 0; i != dimension; ++i) {
        std::string dofname = std::string("Mx") + std::to_string(i);
        doflist.emplace_back(v->toVertex().getDof(dofname));
      }
  return doflist;
} /*}}}*/

}  // namespace geometry

}  // namespace geoxfem
#endif /* DISPLACEMENT_H */
